#include<iostream>
#include<fstream>
#include"PMesh.h"
#include"CSolver.h"

using namespace std;

#define BETA 0.5
#define DT 0.1
#define CP 1950
#define KX 1.73
#define KY 1.73
#define RHO 2520
#define Q 1.2
#define H 1.2
#define W 10
#define GEN 20

int main()
{

	PMesh pm=PMesh(15,12);
	//pm.setNumOfN(121);
	//pm.setNumOfE(202);
	pm.readNodeData("node70x10.txt");
	pm.readElementData("element70x10.txt");
	
	Vector initial=Vector::zeros(pm.NumOfN());
	for(int i=1;i<=pm.Hsize();i++){
		initial(i)=27.0;
	}
	for(int i=1;i<=pm.Vsize()-1;i++){
		initial(i*pm.Hsize()+1)=27.0+5*i;
	}
	for(int i=2;i<=pm.Vsize();i++){
		initial(i*pm.Hsize())=27.0+5*(i-1);
	}
	
	
	Matrix g;
	Vector gv,res;
	pm.setTimeCN(0.5,3600);
	CSolver solv;
	pm.assignInitial(initial);
	ofstream out("bg_temp70x10.dat");
	cout<<"iterasi....\n";
	for(int t=1;t<=240;t++){
		g=pm.assembleStiff(RHO,CP,KX,KY,H,W);
		gv=pm.assembleForce(Q,H,W,GEN,100.0,RHO,CP,KX,KY);
		res=solv.solveGauss(g,gv);
		pm.updateT(res);
		if(t%24==0){
			cout<<"iterasi t= "<<t*3600/(24*3600)<<" days"<<endl;
			out<<"iterasi t= "<<t*3600/(24*3600)<<" days"<<endl;
			for(int i=1;i<=pm.NumOfN();i++){
				out<<pm[i].X()<<"\t"<<pm[i].Y()<<"\t"<<pm[i].T()<<"\n";
			}
		}
	}
	out.close();
	
	return 0;
}

